library(Seurat)
library(scRepertoire)#
library(ggplot2)
library(cowplot)
library(scater)
library(scran)
library(dplyr)
library(Matrix)
library(muscat)#
library(reshape2)
library(celldex)
library(BiocParallel)
library(BiocNeighbors)
library(data.table)
library(DEsingle)
library(stringr)
library(sva)
library(readxl)
library(DESeq2)
# library(DESeq)#
library(pamr)
library(ggpubr)
library(ggraph)
library(gplots)
# library(pca3d)#
library(rgl)
library(scatterplot3d)
library(FactoMineR)
library(ggfortify)
library(useful)
library(tidyverse)
library(kableExtra)
library(xfun)
library(psych)
library(limma)
library(calibrate)
library(pheatmap)
library(ggraph)
library(circlize)
library(scales)
# install.packages("pca3d")
# BiocManager::install("muscat")


setwd("C:/Users/ZFB/Desktop/单细胞生信/GSE212966")
Only_T0 <- readRDS("./data/temp/T_cluster_id_test_1.rds")
Only_T <- Only_T0  
levels(x = Only_T)
sub<- c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", 
        "NKT", "CD8_Tc", "CD8_Te", "CD8_Tem", "CD8_Tex", "CD8_Trm")
Only_T <- subset(Only_T,idents= sub )
levels(x = Only_T)
# [1] "CD4_Tem"  "CD4_Th"   "CD4_Tn"   "CD4_Treg" "CD8_Tc"   "CD8_Te"   "CD8_Tem"  "CD8_Trm"  "NKT"  
table(Idents(Only_T))
DimPlot(Only_T, reduction = "umap", label = TRUE, pt.size = 0.5) 
#细胞及细胞中基因与RNA数量
slotNames(Only_T)
Only_T@assays
dim(Only_T@meta.data)
View(Only_T@meta.data)
#1. T差异分析####
##T分群绘图####
hms_cluster_id<-Only_T
# hms_cluster_id<-readRDS("./data/temp/T_cluster_id_test_1.rds")

sub <- c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", 
         "NKT", "CD8_Tc", "CD8_Te", "CD8_Tem", "CD8_Tex", "CD8_Trm")
dplist <- list()
for (sub_ in sub){
  # sub[1]
  # sub_=sub[1]
  sub_T <- subset(hms_cluster_id, idents = sub_) #细胞数据
  dplist[[sub_]] <- DimPlot(sub_T, reduction = "umap")
  # 保存分群到temp文件夹
  saveRDS(sub_T, file = paste0("./data/temp/", sub_, ".rds"))
  dplist[[sub_]] <- DimPlot(sub_T, reduction = "umap", label = F, pt.size = 0.5)
}
# 绘图
pdf("./data/output/T_sub_批量UMAP图.pdf",width = 12,height = 9)
plot_grid(dplist[['CD4_Tem']],
          dplist[['CD4_Th']],
          dplist[['CD4_Tn']],
          dplist[['CD4_Treg']],
          dplist[['CD8_Tc']],
          dplist[['CD8_Te']],
          dplist[['CD8_Tem']],
          dplist[['NKT']],
          dplist[['CD8_Trm']])
dev.off()



##Only_T GSEA_补_DEG_20240328 ####

Only_T<-readRDS('./data/temp/Only_T.rds')
# a<-Only_T@meta.data
# write.table(a,"./data/temp/a.csv",sep=",")

class(Only_T)
Only_T.sec<-as.SingleCellExperiment(Only_T)
Only_T.sec

table(Only_T@meta.data$tech)
# Normal  Tumor 
# 7984  10331 
group<-factor(c(rep(1,7984),rep(2,10331)))
table(group)

rds <- Only_T
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/Only_T_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/Only_T_DEG_results_classified.csv",sep = ",")


##CD4_Tn_0####
CD4_Tn<-readRDS("./data/temp/CD4_Tn.rds")
a<-CD4_Tn@meta.data
write.table(a,"./data/temp/a.csv",sep=",")

class(CD4_Tn)
CD4_Tn.sec<-as.SingleCellExperiment(CD4_Tn)

table(CD4_Tn@meta.data$tech)
# Normal  Tumor 
# 1373   1669 
group<-factor(c(rep(1,1373),rep(2,1669)))
table(group)

rds <- CD4_Tn
rds<-readRDS('./data/temp/CD4_Tn.rds')

counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  3042
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD4_Tn_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD4_Tn_DEG_results_classified.csv",sep = ",")


##CD4_Tem_0####
CD4_Tem<-readRDS("./data/temp/CD4_Tem.rds")
write.table(CD4_Tem@meta.data,"./data/temp/CD4_Tem_meta.data.csv",sep=",")
CD4_Tem.sec<-as.SingleCellExperiment(CD4_Tem)
CD4_Tem.sec

rds <- CD4_Tem
table(rds@meta.data$tech)
# Normal  Tumor 
# 541   1355  
group<-factor(c(rep(1,541),rep(2,1355)))
table(group)
# 1    2 
# 541 1335 
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  1896
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD4_Tem_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD4_Tem_DEG_results_classified.csv",sep = ",")

##CD4_Treg_0####
CD4_Treg<-readRDS("./data/temp/CD4_Treg.rds")
write.table(CD4_Treg@meta.data,"./data/temp/CD4_Treg_meta.data.csv",sep=",")
CD4_Treg.sec<-as.SingleCellExperiment(CD4_Treg)
CD4_Treg.sec

rds <- CD4_Treg
table(rds@meta.data$tech)
# Normal  Tumor 
# 421   1007 
group<-factor(c(rep(1,421),rep(2,1007)))
table(group)
# 1    2  
# 676   1204 
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  1428
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD4_Treg_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD4_Treg_DEG_results_classified.csv",sep = ",")


##CD8_Te_0####
CD8_Te<-readRDS("./data/temp/CD8_Te.rds")
write.table(CD8_Te@meta.data,"./data/temp/CD8_Te_meta.data.csv",sep=",")
CD8_Te.sec<-as.SingleCellExperiment(CD8_Te)
CD8_Te.sec

rds <- CD8_Te
table(rds@meta.data$tech)
# Normal  Tumor 
# 405   709 
group<-factor(c(rep(1,405),rep(2,709)))
table(group)
# 1    2  
# 676   1204 
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  1114
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD8_Te_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD8_Te_DEG_results_classified.csv",sep = ",")


##CD8_Trm_0####
CD8_Trm<-readRDS("./data/temp/CD8_Trm.rds")
write.table(CD8_Trm@meta.data,"./data/temp/CD8_Trm_meta.data.csv",sep=",")
CD8_Trm.sec<-as.SingleCellExperiment(CD8_Trm)
CD8_Trm.sec

rds <- CD8_Trm
table(rds@meta.data$tech)
# Normal  Tumor 
# 647   553  
group<-factor(c(rep(1,647),rep(2,553)))
table(group)
# 1    2  
# 676   1204 
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  1880
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD8_Trm_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD8_Trm_DEG_results_classified.csv",sep = ",")


##CD8_Tem_0####
CD8_Tem<-readRDS("./data/temp/CD8_Tem.rds")
write.table(CD8_Tem@meta.data,"./data/temp/CD8_Tem_meta.data.csv",sep=",")
CD8_Tem.sec<-as.SingleCellExperiment(CD8_Tem)
CD8_Tem.sec

rds <- CD8_Tem
table(rds@meta.data$tech)
# Normal  Tumor 
# 668   967  
group<-factor(c(rep(1,668),rep(2,967)))
table(group)
# 1    2  
# 668   967  
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  1653
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD8_Tem_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD8_Tem_DEG_results_classified.csv",sep = ",")




##NKT_1#### 
NKT<-readRDS('./data/temp/NKT.rds')
table(NKT@meta.data$tech)
# Normal  Tumor 
# 1218    957 
class(NKT)
group<-factor(c(rep(1,1218),rep(2,957)))

rds<-readRDS('./data/temp/NKT.rds')
counts<-as.matrix(rds@assays$RNA@counts)
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/NKT_DEG_results.csv",sep = ',')
write.table(results.classified,"./data/temp/NKT_DEG_results_classified.csv",sep = ',')


##CD4_Th_1####
CD4_Th<-readRDS("./data/temp/CD4_Th.rds")
write.table(CD4_Th@meta.data,"./data/temp/CD4_Th_meta.data.csv",sep=",")
CD4_Th.sec<-as.SingleCellExperiment(CD4_Th)
CD4_Th.sec

rds <- CD4_Th
table(rds@meta.data$tech)
# Normal  Tumor 
# 676   1204  
group<-factor(c(rep(1,676),rep(2,1204)))
table(group)
# 1    2  
# 676   1204 
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  1880
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD4_Th_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD4_Th_DEG_results_classified.csv",sep = ",")


##CD8_Tc_1####
CD8_Tc<-readRDS("./data/temp/CD8_Tc.rds")
write.table(CD8_Tc@meta.data,"./data/temp/CD8_Tc_meta.data.csv",sep=",")
CD8_Tc.sec<-as.SingleCellExperiment(CD8_Tc)
CD8_Tc.sec

rds <- CD8_Tc
table(rds@meta.data$tech)
# Normal  Tumor 
# 676   1204  
group<-factor(c(rep(1,676),rep(2,1204)))
table(group)
# 1    2  
# 676   1204 
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  1880
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD8_Tc_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD8_Tc_DEG_results_classified.csv",sep = ",")

##CD8_Tc_new####
# 20240201
# CD8_Tc -> CD8_Tc+CD8_Te
# CD8_Te -> CD8_Tex
CD8_Tc<-readRDS("./data/temp/CD8_Tc.rds")
write.table(CD8_Tc@meta.data,"./data/temp/CD8_Tc_meta.data.csv",sep=",")
CD8_Tc.sec<-as.SingleCellExperiment(CD8_Tc)
CD8_Tc.sec

rds <- CD8_Tc
table(rds@meta.data$tech)
# Normal  Tumor 
# 1361    772  
group<-factor(c(rep(1,1361),rep(2,772)))
table(group)
# 1    2  
# 676   1204 
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  2133
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD8_Tc_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD8_Tc_DEG_results_classified.csv",sep = ",")


##CD8_Te_new####
# 20240201
# CD8_Tc -> CD8_Tc+CD8_Te
# CD8_Te -> CD8_Tex
CD8_Te<-readRDS("./data/temp/CD8_Te.rds")
write.table(CD8_Te@meta.data,"./data/temp/CD8_Te_meta.data.csv",sep=",")
CD8_Te.sec<-as.SingleCellExperiment(CD8_Te)
CD8_Te.sec

rds <- CD8_Te
table(rds@meta.data$tech)
# Normal  Tumor 
# 674   1138  
group<-factor(c(rep(1,674),rep(2,1138)))
table(group)
# group
# 1    2 
# 674 1138 
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  1812
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD8_Te_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD8_Te_DEG_results_classified.csv",sep = ",")

##CD8_Tex_new####
# 20240201
# CD8_Tc -> CD8_Tc+CD8_Tex
# CD8_Te -> CD8_Tex
CD8_Tex<-readRDS("./data/temp/CD8_Tex.rds")
write.table(CD8_Tex@meta.data,"./data/temp/CD8_Tex_meta.data.csv",sep=",")
CD8_Tex.sec<-as.SingleCellExperiment(CD8_Tex)
CD8_Tex.sec

rds <- CD8_Tex
table(rds@meta.data$tech)
# Normal  Tumor 
# 674   1138  
group<-factor(c(rep(1,405),rep(2,709)))
table(group)
# group
# 1    2 
# 405 709
counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  1114
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,"./data/temp/CD8_Tex_DEG_results.csv",sep = ",")
write.table(results.classified,"./data/temp/CD8_Tex_DEG_results_classified.csv",sep = ",")




##volcano_plot####
###法1🪲####
res <- read.csv("./data/temp/CD4_Tn_DEG_results.csv", header=TRUE,sep=",")
head(res)

library(tidyverse)
library(ggrepel)
res$logFC <- -log2(res$norm_foldChange)#⚠️负号是因为之前Normal/Tumor
res$P.Value <- res$pvalue
res$adj.P.Val <- res$pvalue.adj.FDR

data <- 
  res %>% 
  mutate(change = as.factor(ifelse(pvalue.adj.FDR < 0.05 & abs(norm_foldChange) > 1,
                                   ifelse(norm_foldChange > 1 ,'Up','Down'),'No change'))) %>% 
  rownames_to_column('gene')
# head(data)


core_gene <- c("PRSS1","CLPS","CTRB1","PNLIP","CPA1","CELA3A","HSP90AB1","PRSS2","TXNIP","CPB1",
               "HSP90AA1","DNAJA1","HSPD1","HSPE1","CTRC","PLA2G1B","TRB2","ELA3B","YCN","AMY2A")
core_gene

ggplot(data,aes(logFC, -log10(adj.P.Val)))+
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black")+
  geom_vline(xintercept = c(-1.2,1.2), linetype = "dashed", color = "black")+
  geom_point(aes(size = -log10(adj.P.Val), 
                 color = -log10(adj.P.Val)))+
  scale_color_gradientn(values = seq(0,1,0.2),
                        colors = c("#39489f","#39bbec","#f9ed36","#f38466","#b81f25"))+
  scale_size_continuous(range = c(0,1))+
  theme_bw(base_size = 15)+
  theme(panel.grid = element_blank(),
        legend.position = 'right',
        legend.justification = c(0,1))+
  # 设置图例
  guides(col = 
           guide_colorbar(title = "-Log10_q-value",
                          ticks.colour = NA,
                          reverse = T,
                          title.vjust = 0.8,
                          barheight = 8,
                          barwidth = 1),
         size = "none") +
  # 添加标签：
  geom_text_repel(data = filter(data, gene %in% core_gene),
                  max.overlaps = getOption("ggrepel.max.overlaps", default = 20),
                  # 这里的filter很关键，筛选你想要标记的基因
                  aes(label = gene),
                  size = 2, 
                  color = 'black') +
  xlab("Log2FC")+
  ylab("-Log10(FDR q-value)")+
  # x y 轴显示范围和刻度间隔
  scale_x_continuous(limits =c(-10,10),breaks =seq(-10,10,2))+
  scale_y_continuous(limits =c(0,20),breaks =seq(0,20,2))

###法2####
res <- read.csv("./data/temp/CD4_Tn_DEG_results.csv", header=TRUE,sep=",")
title <- paste0("CD4_Tn"," Volcano Plot")
res$log2FoldChange <- -log2(res$norm_foldChange)
#⚠️这里取-负数是因为Normal/Tumor比反了
res$fdr <- res$pvalue.adj.FDR

a1 <- subset(res, fdr<.05 & abs(log2FoldChange)>1)
a2 <- subset(res, fdr<.05 & log2FoldChange>1)

# head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main=title, xlim=c(-6,6),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
# with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
# with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

res <- read.csv("deg_Effector_Memory_CD4_T", header=TRUE,sep="\t")
head(res)
with(res, plot(log2FoldChange, -log10(fdr), pch=20, main="Volcano plot", xlim=c(-6,7),col="grey"))
# Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
#with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
#with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
abline(h=1.3,v=1,lty=3)
abline(v=-1,lty=3)

##循环出图####
# sub = c("CD4_Tem", "CD4_Th", "CD4_Tn", "CD4_Treg", "CD8_Tc",
#         "CD8_Te", "CD8_Tem", "NKT","CD8_Trm" )# 修改

sub = c("CD8_Tc", "CD8_Te", "CD8_Tex", "NKT")# 修改

for(sub_ in sub){
  # sub_ <- sub[1]
  read_path <- paste0("./data/temp/",sub_,"_DEG_results.csv")
  save_path <- paste0("./data/output/",sub_,"_VolcannoPlot.pdf")
  read_path
  res <- read.csv(read_path, header=TRUE,sep=",")
  title <- paste0(sub_," Volcano Plot")
  #⚠️⚠️⚠️⚠️⚠️⚠️⚠️⚠️
  res$log2FoldChange <- -log2(res$norm_foldChange)  #⚠️这里取-负数是因为Normal/Tumor
  #⚠️⚠️⚠️⚠️⚠️⚠️⚠️⚠️
  res$fdr <- res$pvalue.adj.FDR
  pdf(save_path, width = 10, height = 10)
  with(res, plot(log2FoldChange, -log10(fdr), pch=20, 
                 main="", xlim=c(-6,6),col="grey"))
  # Add colored points: red if pvalue<0.05, orange of log2FC>1, green if both)
  # with(subset(res, fdr<.05 ), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
  # with(subset(res, abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="orange"))
  with(subset(res, fdr<.05 & abs(log2FoldChange)>1), points(log2FoldChange, -log10(fdr), pch=20, col="blue"))
  with(subset(res, fdr<.05 & log2FoldChange>1), points(log2FoldChange, -log10(fdr), pch=20, col="red"))
  #with(subset(res, P.Value<.05 & abs(log2FC)>1), textxy(logFC, -log10(P.Value), labs=Gene, cex=.6))
  abline(h=1.3,v=1,lty=3)
  abline(v=-1,lty=3)
  
  title(main = list(title, cex = 3, col = "black"))
  dev.off()
}


##to_KMplot_获取筛选基因列表，用于####
#这里只筛选了上调的部分基因，关注上调基因取KMplot寻找有生存差异的基因
sub = c("CD4_Treg", "CD8_Tex", "NKT")# 修改
for(sub_ in sub){
  # sub_ <- sub[1]
  read_path <- paste0("./data/temp/",sub_,"_DEG_results.csv")
  save_path <- paste0("./data/output/",sub_,"_DEG_results_to_KMplot.csv")
  read_path
  save_path
  res <- read.csv(read_path, header=TRUE,sep=",")
  
  #⚠️⚠️⚠️⚠️⚠️⚠️⚠️⚠️
  res$log2FoldChange <- -log2(res$norm_foldChange)  #⚠️这里取-负数是因为Normal/Tumor
  #⚠️⚠️⚠️⚠️⚠️⚠️⚠️⚠️
  res$fdr <- res$pvalue.adj.FDR
  #pvalue和foldchange都满足的
  res <- subset(res, fdr<.05 & log2FoldChange>1)#这里只筛选了上调的部分基因
  res <- res[,c("norm_foldChange","log2FoldChange","pvalue.adj.FDR")]
  write.table(res,save_path,sep=',')
  
}
## to GSEA获取 20240327 ####
###1.to_GSEA Tsub ####
#单独细胞亚群GSEA分析
#重新分群之后，一部分内容放在data/temp2，现在全放在了data/temp
#classified 里面有上调的，在toKMplot其实也可以选择里面的down，GSEA这里上下调都保留，不做特别的筛选
#补充：GSEA分析数据不用差异显著的数据
sub=c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", 
      "NKT", "CD8_Tc", "CD8_Te", "CD8_Tem", "CD8_Tex", "CD8_Trm")
for(sub_ in sub){
  # sub_ <- sub[1]
  read_path <- paste0("./data/temp/",sub_,"_DEG_results_classified.csv")
  save_path <- paste0("./data/temp/",sub_,"_DEG_results_to_GSEA.csv")
  read_path
  save_path
  res <- read.csv(read_path, header=TRUE,sep=",")
  
  #⚠️⚠️⚠️⚠️⚠️⚠️⚠️⚠️
  res$norm_foldChange <- 1/res$norm_foldChange
  res$log2FoldChange <- log2(res$norm_foldChange)
  # res$log2FoldChange <- -log2(res$norm_foldChange)  #⚠️这里取-负数是因为Normal/Tumor
  #⚠️⚠️⚠️⚠️⚠️⚠️⚠️⚠️
  res$fdr <- res$pvalue.adj.FDR
  #pvalue和foldchange都满足的
  # res <- subset(res, fdr<.05 & (log2FoldChange>1 | log2FoldChange< -1))
  res <- res[,c("norm_foldChange","log2FoldChange","pvalue.adj.FDR")]
  write.table(res,save_path,sep=',')
  
}
### 2.to_GSEA Only_T####
sub_ <- "Only_T"
read_path <- paste0("./data/temp/",sub_,"_DEG_results_classified.csv")
save_path <- paste0("./data/temp/",sub_,"_DEG_results_to_GSEA.csv")
read_path
save_path
res <- read.csv(read_path, header=TRUE,sep=",")

#⚠️⚠️⚠️⚠️⚠️⚠️⚠️⚠️
res$norm_foldChange <- 1/res$norm_foldChange
res$log2FoldChange <- log2(res$norm_foldChange)
# res$log2FoldChange <- -log2(res$norm_foldChange)  #⚠️这里取-负数是因为Normal/Tumor
#⚠️⚠️⚠️⚠️⚠️⚠️⚠️⚠️
res$fdr <- res$pvalue.adj.FDR
#pvalue和foldchange都满足的
# res <- subset(res, fdr<.05 & (log2FoldChange>1 | log2FoldChange< -1))
res <- res[,c("norm_foldChange","log2FoldChange","pvalue.adj.FDR")]
write.table(res,save_path,sep=',')




#2. Fibroblast差异分析####
##2.1 iCAF####
# 有时间再跑把，太费时间了
pbmc <- readRDS(sub_cluster_id, file = "./data/temp/Fibroblast_3亚群_i-m-myCAF.rds")
sub_="iCAF"
sub <- subset(pbmc, idents = sub_) #细胞数据
paste0("./data/temp/",sub_,"_meta.data.csv")
write.table(sub@meta.data,paste0("./data/temp/",sub_,"_meta.data.csv"),sep=",")
sub.sec<-as.SingleCellExperiment(sub)
sub.sec

rds <- sub
table(rds@meta.data$tech)
# Normal  Tumor 
# 1336   3244  
group<-factor(c(rep(1,1336),rep(2,3244)))
table(group)

counts<-as.matrix(rds@assays$RNA@counts)
dim(counts)
#  27579  4580
results<-DEsingle(counts=counts,group=group)
results.classified <- DEtype(results = results, threshold = 0.05)
write.table(results,paste0("./data/temp/",sub_,"_DEG_results.csv"),sep = ",")
write.table(results.classified,paste0("./data/temp/",sub_,"_DEG_results_classified.csv"),sep = ",")

